from brian import *
from brian.library.random_processes import *
from brian.library.synapses import  *
import mydataRG as mydatrg

md = mydatrg.RGVars()
md.__init__()
#defaultclock.dt = 0.005*ms

def i_inj(t):
    if t <= 1.1*ms and t > 1 *ms:
        return 0*nA
    else:
        return 0*nA
    
    
DUR = 200

def spToList(a):
    l = list()
    
    for x in a.spiketimes:
        l+=((a[x]/20*1000/md.BINSIZE)/ms).tolist()
    return l

def napvtrap(x,y):
    a=x/y
    if abs(a)<0.0001:
        return  1/(2+a)
    else:
        return 1/(exp(a)+1)
    #return 1/(exp(a)+1)
    
def napefun(x,y):
    a = x/y 
    if abs(a)<0.0001:
        return 1+a
    else:
        return exp(a)
    #return exp(a)
            
 
def m_inf_na(V):
    return napvtrap(-(V+35),7.8)   
    
def h_inf_na(V):
    return napvtrap((V+55),7)

def tau_h_na(V):
    return msecond*(30/(napefun((V+50),15)+napefun(-(V+50),16)))


def tau_m_k(V):
    return msecond*(7/(napefun(V+40,40)+napefun(-(V+40),50)))

def m_inf_nap(V):
    return napvtrap(-V-47.1,3.1)

def h_inf_nap(V):
    return napvtrap(V+59,8)

def tau_h_nap(V):
    return msecond*(2*1200/(napefun(V+59,16)+napefun(-V-59,16)))


def current_inj(t):
    if t > 0*msecond:
        return -0.5*namp
    else:
        return 0*mamp

# The model
eqsPF=Equations('''
dvm/dt = (i_inj(t)-(md.gna*mna*mna*mna*hna*(vm-md.Ena))-(md.gk*mk*mk*mk*mk*(vm-md.Ek))-(md.gnap*mnap*hnap*(vm-md.Ena))-md.gl*(vm-md.El))/md.cpf : mvolt
dmna/dt = (m_inf_na(vm/mvolt)-mna)/md.taumna : 1
dhna/dt = (h_inf_na(vm/mvolt)-hna)/tau_h_na(vm/mvolt) : 1
dmk/dt = (m_inf_k(vm/mvolt)-mk)/tau_m_k(vm/mvolt) : 1
dmnap/dt = (m_inf_nap(vm/mvolt)-mnap)/md.tau_m_nap : 1
dhnap/dt = (h_inf_nap(vm/mvolt)-hnap)/tau_h_nap(vm/mvolt) : 1

''')


P=NeuronGroup(1,model=eqsPF,
threshold=EmpiricalThreshold(threshold=-40*mV),
implicit=True)

# Initialization

vv =md.vinit
P.vm=vv*mV 
P.mna = m_inf_na(vv)
P.hna = h_inf_na(vv)
P.mk = m_inf_k(vv)

P.mnap = m_inf_nap(vv)
P.hnap = h_inf_nap(vv)
#P.El = -64*mV * randn(len(P))*0.64*mV


tracePF=StateMonitor(P,'vm',record=[0])

run(DUR*msecond)
figure(1)

plot(tracePF.times/ms,tracePF[0]/mV)

show()